Theory of antibound states in partially filled narrow band systems 
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We present a theory of the dynamical two-particle response function in the Hubbard model based 
on the time-dependent Gutzwiller approximation. The results are in excellent agreement with 
exact diagonalization on small clusters and give reliable results even for high densities, where the 
usual ladder approximation breaks down. We apply the theory to the computation of antibound 
states relevant for Auger spectroscopy and cold atom physics. A special bonus of the theory is its 
computational simplicity. 
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Much of our understanding of strongly correlated elec- 
tronic systems comes from dynamical responses like the 
one-particle spectral function which, under certain ap- 
proximations, is probed by photoemission and inverse 
photoemission experiments. Less explored is the two- 
particle spectral function in which one studies how the 
system responds to the addition or removal of two par- 
ticles. In the case of two holes in an otherwise filled 
band, the response was computed exactly by Cini[l|, Q 
and SawatzkyQ (CS) in connection with Auger spec- 
troscopy. They showed that for strong enough on-site 
repulsion the spectral function gets dominated by anti- 
bound states in which the two holes propagate together 
paying a large Coulomb cost. 

Despite the interest of the problem, the dynamical two- 
particle response and the formation of antibound states 
in partially-filled correlated systems are not well under- 
stood. Cini and collaborators [J @| have compared ap- 
proximations for the spectral function developed by sev- 
eral groups with exact diagonalization on finite clusters. 
They observed that any attempt to improve the single- 
fcrmion propagators with self-energy corrections or mak- 
ing them self consistent leads to worse results due to the 
lack of vertex corrections which, if included, would tend 
to "undress" the Green's functions. Thus for small filling, 
the best approximation corresponds to a trivial general- 
ization of the original theory, namely summing a ladder 
series with bare Green's functions. For moderate filling 
and for large interactions, this bare ladder approximation 
(BLA) breaks down and no reliable theory is available. 

Several effects are expected to be relevant in the case 
of a partially-filled band. First, strong correlation pro- 
duces band narrowing, which should help to split-off an- 
tibound states from the two-particle continuum. Second, 
the spectral weight of the antibound state should depend 
on doping, since the probability to find an empty site 
where to create an antibound pair depends on the fill- 
ing. Third, the other holes present in the system are 
expected to screen the effective interaction among the 



added holes, which may lead to a renormalization of the 
position of the antibound state with respect to the con- 
tinuum. Last, the chemical potential has a jump as a 
function of doping across the Mott insulating phase of 
narrow-band systems, which should show up in the po- 
sition of the two-particle continuum with respect to the 
antibound state. 

In this work, we present a theory of antibound states 
for the Hubbard model, which incorporates these ef- 
fects. It is based on the computation of pairing fluc- 
tuations within the time-dependent Gutzwiller approx- 
imation (TDGA)0j. Our approach reproduces the ef- 
fects discussed above, while keeping the simplicity of CS 
theory. Interestingly, we find that the effect of a finite 
density is to antiscreen the Hubbard U interaction, i.e., 
the effective interaction is larger than the bare one and 
becomes singular as the Mott phase is approached [c.f. 
Fig. QJa)] . The comparison of our results with exact di- 
agonalizations shows that TDGA is reliable even at high 
densities where the BLA breaks down (c.f. Fig. [3]). 

Our starting point is the Hubbard Hamiltonian: 
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where cj creates a fermion with spin a at site i, n{ a = 
c\ a Ci a , U is the on-site repulsion, denotes the hopping 
amplitude, and /i is the chemical potential. 

We arc interested in the following two-particle response 
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for u; > (lu < 0) the imaginary part of Eq. $2$ gives the 
two-particle addition (removal) spectra. For the Auger 
application one should consider other effects which have 
been extensively discussed in the literature [1] and will not 
be treated here (e.g., finite life time of the core holeQ 
and the interaction of the core hole with the valence 
electrons 0). 
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FIG. 1: Panels (a), (c), and (d) are obtained within a 
model with a flat density of states and a bare bandwidth 
W {Ubr = 2W) in HF (dashed line) and GA (full lines), (a) 
Effective particle-particle interaction V for different dopings 
n — 1. Negative and positive dopings coincide, (b) Sketch of 
the principal energy scales of lmPu. (c) Energy distance be- 
tween the center of the two-particle scattering states (2E) and 
the doublon energy U for different fillings. The 1 + , 1~ fillings 
are infinitesimal deviations from half filling and coincide in HF 
and in GA for U < Ubr- The same plot represents the effec- 
tive interaction rescaled by the doping V(n — 1) [see Eq. ©]. 
(d) Boundaries of the two-particle continuum uj' = 2S ±ZqW 
for U = 2U br as a function of band filling. 
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Our approach is 
function@,[[o|: |$) = P g \4>) where P, 
out doubly occupied sites from 
to be a Bogoliubov vacuum. We define the single- 
particle density matrix pi a ja> = {4>Wja lCi ° W) an( ^ P a ^ r 
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ing constraints [llj 
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The first step is to construct the charge rotationally in- 
variant energy functional E = in the Gutzwiller 
approximation (GA). This is more easily done by rotating 
at each site the fermion annihilation and creation oper- 
ators to a basis where the anomalous expectation values 
vanish[12]. Then, one derives the GA with one of the 
known techniques 13, LL4[ and rotates back to the origi- 
nal operators. Restricting to a paramagnetic state one 
finds 

E\p, K, D] = ^ tijZiZjPivja + U^Di, (4) 
ija i 

with the hopping renormalization factors 



Di + Ji z (V Di — Ji z — Ji + y/Di~— Ji z + Ji 



Here we defined J ix = (k^ai + k* t ^)/2, J iy 
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double occupancy Di — ($|n,-f7i,-J$). The ground state 
is found by minimizing Eq. ^ with the constraints ([3]), 
leading to the static p°, K°, J° and D°. We will consider 
a paramagnetic normal metal thus n° = J° = J® = 0. 

To compute the response function we add a weak time- 
dependent pairing field F(t) = ~}^i{fi e ~ lultc il c i]+h.c.) to 
Eq. (UJ). This produces small time-dependent deviations 
Sp(t) = p(t) — p°. In addition, since F does not conserve 
the particle number, it induces pairing correlations k, 
which we compute in linear response. 

Previously [E El El, 

the energy was expanded to sec- 
ond order in terms of particle-hole fluctuations, leading 
to effective matrix elements for charge and spin excita- 
tions. For a normal paramagnet neither those channels 
nor SD fluctuations mix with the particle-particle chan- 
nel which simplifies the formalism. The remaining part 
follows text book computations in nuclear physics [llj. 

Expanding the energy up to second order in 8p and k 
one finds: 

^ = ^-^ + 7^(4 + 4). (5) 
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Here £k = z^e^ + Eg denotes the GA dispersion relation 
(ek is the bare one), Eg coincides with the Lagrange 
parameter of the slave boson method[l3| and is given by 
Eg = z z' e with e = J^ia^jPiaja^ z o is the hopping 
renormalization factor at the saddle point and z' is its 
density derivative. Our notation emphasizes the fact that 
Eg can be interpreted as a local GA self-energy. Finally, 
the effective on-site particle-particle interaction is 
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where n denotes the particle concentration. At half fill- 
ing (n = 1), both the numerator and the denominator 
tend to zero and one finds V — U (1 — U/2Ubr) (1 + 
U/Urr)/(\—U/Urr), which coincides with the particle- 
hole casejy, [lfj, [l5|. Here Ubr = 8e is the criti- 
cal interaction for the Brinkman-Rice metal insulator 
transition [Tol El] . 

The response function can be readily derived from the 
equations of motion of the pair matrix in a normal system 
after using the constraint (|3|) to express the first term in 
Eq. ([5]) as a quadratic contribution in k The momen- 
tum dependent pair-correlation function is given by the 
usual ladder expression but with the effective interaction 
of Eq. ©: 
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where P (q, ui) is the non-interacting two-quasiparticle 
correlation function 
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FIG. 2: Local spectral function for different fillings in TDGA 
and BLA. Results are for the Hubbard model on a square lat- 
tice with nearest-neighbor hopping (Ubr = 128t/-7r 2 ). The 
vertical arrows indicate the position of 2/j,, separating the ad- 
dition part &' > 2[i from the removal part uj' < 2fi. The 
intensity of two-particle scattering states have been multi- 
plied by 10 3 . Inset: the n = case which coincides in the two 
approximations. 



evaluated with the GA dispersion relation £k- N s is the 
number of sites, f(£u.) is the Fermi distribution function, 
and ?7k,k' = 0+sign(e k + £k' - 2/i). 

Eqs. ([6]), ([7]), and ([8]) constitute our main result. Our 
approach leads to the same formal ladder structure as 
in the CS theory[H, H E3] but with the HF self-energy 
(T,hf = Un/2, zq = 1) replaced by the GA one and the 
Hubbard repulsion U replaced by an effective interaction 
V. Notice that the "new" Eq. © is valid in the BLA 
provided one replaces Eg — > ?->hf leading to V = U. 

Fig. QJa) shows V [c.f. Eq. ((6|)] as a function of U for 
a band of bare width W and a flat density of states. For 
fillings close to or 2 or small U, the effective interac- 
tion is close to the bare U, as expected. By contrast, as 
the Mott phase is approached, V diverges. This singular 
behavior is essential for the correct description of dense 
systems close to the Mott insulator. 

For the following analysis, it is convenient to shift the 
origin of energies to eliminate the chemical potential in 



Eq. ©, i.e., we define u' = uj + 2/i. Following Refs. [18|, 
Il9j . one can compute exactly P(Q,w) at Q = (tt, ..,7t) 
for the Hubbard model with nearest-neighbor hopping. 
In this case, the full spectral function is exhausted by a 
single pole at uj' = U. The antibound state consist of a 
doublon, i.e., an on-site pair. This provides a quick and 
instructive check of the theory. Indeed, by using that 
ek + ek+Q = 0, one can verify that both BLA and the 
present theory reproduce the exact result. For general 
momenta and large U, Eq. ([7]) has a single pole for uj' ^ U 
(i.e., the antibound state) and a continuum at low (high) 
energy for n < 1 (n > 1). The local response is obtained 
as Pu(uj) — 1/N S -P(q, uj). Fig. QJb) shows a sketch 
of ImPa(uj') with the continuum at 2£ — ZqW < uj' < 



2£ + ZqW and the antibound state at uj' ~ JJ. The dip 
in the continuum at 2\x separates the addition part for 
uj' > 2/i from the removal part for uj' < 2/i. 

The problem within the BLA is not so much the en- 
ergy of the antibound state but rather the position of 
the two-particle continuum, which is given by the HF 
eigenvalues. This affects the antibound state because, as 
the continuum approaches the energy U , the antibound 
pair becomes less localized in the relative coordinate and 
eventually disappears. The distance between the contin- 
uum (~ 2S) and the antibound state (~ U) as a function 
of U is shown in Fig.QJc). The picture can be easily un- 
derstood by noticing that by rescaling the y axis by 1/2 
one obtains one-particle energies. For an almost filled or 
an almost empty band, as well as for small U, the posi- 
tion of the HF and the GA band coincide. At half filling, 
the HF self-energy is T, HF = U/2 so that 2S - U = 0. 
This coincides with the GA for U < Ubr, however, for 
U > Ubr the GA self-energy bifurcates in two solutions, 
corresponding to infinitesimal positive and negative de- 
viation from half filling, due to the opening of the Mott- 
Hubbard gap. Thus for U larger than U br and moderate 
filling, the HF band is close to U, whereas the GA band 
is well separated from it. In this case we can anticipate 
quite different two-particle spectra in the two approxi- 
mations. This dramatic difference is also illustrated in 
Fig. [ljd) where the boundaries of the continuum with 
respect to the doublon energy are shown as a function of 
filling for U = 2Ub r- The HF self-energy leads to a linear 
evolution. By contrast, in GA the band remains nearly at 
the same energy and narrows when n — > 1 due to correla- 
tion. At n = 1 the band jumps due to the Mott-Hubbard 
gap and the situation reverses. Clearly the GA contin- 
uum is always far from the antibound state whereas, in 
HF, it is generally much closer and overlaps the uj' = U 
line in a large range of filling near n = 1. Therefore, the 
formation of tight antibound states will be much more 
favored in the GA case. 

Fig. [2] compares the local two-particle spectral function 
for an infinite two-dimensional system and n < 1, within 
TDGA and BLA. The inset shows the n = case where 
TDGA and BLA coincide. Differences occur at finite 
concentrations (main panel) where the line shapes are 
dominated by the antibound state at uj' ~ U (as in CS), 
which is significantly stronger in the TDGA. The inten- 
sity of the continuum at low energies has been multiplied 
by 10 3 to make the line shape visible. As anticipated the 
two-particle continuum is far from the antibound state 
in GA, whereas it quickly approaches it in the BLA. 

The antibound state can propagate and forms a band 
which gives the width of the high-energy feature. The 
lower edge of this band corresponds to q = (tt, tt) for 
n < 1 and is at uj' = U. For large U, the bandwidth 
is of order t 2 /U for n = 0[3] but becomes of order t 
(specifically 2zg|e/(l — n)\) for finite n, since the kinetic 
energy can move a doublon at first order if there is a 
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single occupied site next to it. 

The pair correlation function satisfies the sum rules: 

1 f°° 

/ dwlmPuicj) = l-n + {n^ni\) (9) 

1 f 2 ^ 

/ du'lm Pii(w') = (ni-i-nii). (10) 

T J-oo 

This can be used to evaluate ladder corrections to the 
GA or HF double occupancy. The fact that the area 
of the removal part is much larger in the BLA than in 
TDGA reflects a larger double occupancy in the former. 
This is not surprising since at zero order BLA neglects 
correlations at all. Furthermore, in TDGA, as the system 
approaches the Mott phase, the hopping renormalizes to 
zero and the system becomes more "atomic" like. This 
explains the vanishing of two-particle scattering states as 
n — > 1. (Clearly in an exact computation, a small finite 
double occupancy and scattering intensity will remain 
in the Mott phase). Contrary, in the BLA the system 
becomes more "band" like as the filling is increased due 
to the closing of the gap between the scattering states and 
the doublon energy. Indeed for n = 0.85 the antibound 
state exists only for some values of the momentum. 

In order to validate our results, we have computed the 
exact two-particle addition spectra for 10 particles on a 
4x4 lattice with only nearest-neighbor hopping t and 
U/t = 15, by using exact diagonalization. Fig. [3]shows a 
comparison between the present theory and BLA. Here, 
we go back to our original variables and fix the origin 
of energy at 2fi. Despite the large value of the Hubbard 
repulsion, TDGA yields excellent agreement with exact 
diagonalization concerning the location, width and inten- 
sity of the high-energy antibound states. On the other 
hand, BLA predicts that these excitations have a much 
lower energy when referenced to 2fj, and no clear separa- 
tion with the band states is visible (see upper- left inset). 
For the system under consideration, there are three band- 
like two-particle energies which are very well reproduced 
by TDGA in contrast with BLA. The upper-right inset 
demonstrates that the double occupancy after Eq. (J9j) is 
accurate within TDGA, whereas BLA overestimates it as 
expected. The excellent performance of TDGA is not re- 
stricted to this particular value of U but persists to even 
larger (and of course lower) on-site repulsions. 

To conclude, we have presented a computation of pair 
fluctuations for the Hubbard model exhibiting antibound 
states for large Coulomb repulsion. Our approximation 
gives reliable results even for large densities, where we 
are not aware of any accurate theory. The simplicity 
of the method suggests its application to the compu- 
tation of Auger spectra on top of realistic Gutzwillcr 
calculations [20j . Our theory can also be applied to ultra- 
cold fcrmion atoms in optical lattices, which can be de- 
scribed by the Hubbard model as well[21]. The possibil- 
ity to observe antibound states has already been demon- 
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FIG. 3: Imaginary part of the pair correlation function for 
the Hubbard model with U/t = 15 and 10 particles on a 
4x4 square lattice obtained by exact diagonalization, TDGA 
and BLA. The origin of energy is at 2/j,. The upper-left in- 
set enlarges the region of low-energy band excitations and 
the upper-left inset shows the frequency evolution of the in- 
tegrated spectra. The broadening of the delta-peaks is 0.5t 
in the main panel and O.lt in the insets. The arrow indicates 
the exact U — 2/i value. 

strated in the Bose case|22j|. 
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